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Nonadiabatic unitary evolution with tailored time-dependent Hamiltonians can prepare systems of cold 
atomic gases with various desired properties. For a system of two one-dimensional quasicondensates cou- 
pled with a time-varying tunneling amplitude, we show that the optimal protocol, for maximizing any figure of 
merit in a given time, is bang-bang, i.e., the coupling alternates between only two values through a sequence of 
sudden quenches. Minimizing the energy of one of the quasicondensates with such nonadiabatic protocols, and 
then decoupling it at the end of the process, can result in effective cooling beyond the current state of the art. 
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Recent advances in the physics of ultracold atoms, brought 
about by ingenious cooling techniques such as evaporative and 
laser cooling, have stirred up great interest in the nonequi- 
librium dynamics of many -body quantum systems uJ-lsl]. 
Despite remarkable progress, however, the quest for lower 
temperatures still continues. For example, creating ground 
states of important model Hamiltonians, such as the two- 
dimensional fermionic Hubbard model, has remained elusive. 
In addition to cooling, preparing systems with other desired 
characteristics, such as, e.g., number-squeezed ones, is of con- 
siderable interest due to the potential applications to quantum 
metrology and precision measurements. Focusing on a a pair 
of two coupled elongated (assumed one-dimensional) quasi- 
condensates (hereafter referred to simply as condensates de- 
spite lack of true long-range order), we propose a scheme for 
preparing cold atomic systems with custom-ordered figures of 
merit through optimal control of their nonequilibrium quan- 
tum dynamics. As we will show, the large degree of dynami- 
cal control over these systems provides, among others, a new 
means of bringing them even closer to zero temperature. 

Let us begin by giving a few examples of experimentally 
relevant quantities one can optimize in cold atom systems: 

1. Eff'ective cooling: by minimizing quantities such as the 
excess energy, number of quasiparticle excitations, or 
the trace distance between the density matrix of the sys- 
tem and its zero-temperature density matrix. 

2. Phase coherence: by minimizing the fluctuations of the 
relative phase between two condensates (with spatially 
fluctuating phases'), which is important for matter-wave 
interferometry ll6t-[l0ll. 

3. Number squeezing: by minimizing the particle number 
fluctuations of a system fll ll Il2 L which is important, 
e.g., in precision measurements ll3 l. 



Focusing on effective cooling, we show in this paper that (at 
least) one of the condensates in the system of Fig. [T]can be 
cooled down by a factor of 5 with our proposed method under 
reasonable experimental conditions. This is not a fundamental 
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FIG. 1: Two coupled one-dimensional condensates. The tunneling 
amplitude A(f) can be tuned by raising or lowering the potential bar- 



bound, however, and cooling by several orders of magnitude 
is in principle possible for highly asymmetric systems. 

Let us now formulate our proposal generically. Consider 
a quantum system with Hamiltonian H, which is comprised 
of two coupled subsystems: H - Hi + H2 + V, where //, is 
the Hamiltonian of subsystem / and V is the coupling Hamil- 
tonian. Generically, V is a sum of certain local terms, with 
some coupling constants {A}. We assume that i) we have 
time-dependent control over the coupling constants {A{t)], i.e., 
within a range determined by the experimental constraints, we 
can tune them to any value as a function of time, ii) for all 
initial [Aq], we can prepare the system at inverse temperature 
y6o with the current state-of-the-art cooling methods, and, iii) 
we have a fixed time t to carry out a dynamical process (by 
tuning the Hamiltonian), during which the system undergoes 
quantum-coherent unitary evolution. 

Our goal is to find an optimal protocol {A{t)} such that, 
at the end of the process (f = r), the energy of subsys- 
tem 1, or some other custom-ordered cost function, is min- 
imized. For a given protocol, the density matrix evolves as 
p(f) — i{H{{A{t)]), pit)] with initial conditions determined by 
the thermal state at f = 0. Thus, p(r), and, consequently, cost 
functions such as £i(t) - tr [H\p{t)\ are functionals of {A{t)], 
< f < T. Notice that if we decouple the two subsystems at 
time T, the energy of subsystem 1 will remain equal to £i(t) 
for all subsequent times. 

The key question addressed in this paper is how to minimize 
this functional of {/1(f))- We find that i) the Pontryagin's max- 
imum principle provides a deep understanding of the struc- 
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ture of such protocols, and ii) the simulated annealing method 
used in Ref. lll4ll gives a simple and generic way for perform- 
ing such optimization. In simulated annealing, we discretize 
time, approximate an arbitrary protocol by a piece-wise con- 
stant one, and perform direct (classical) Monte-Carlo (MC) 
simulations with kinetic moves consisting of small random 
displacements of randomly chosen pieces of the protocol ll4ll . 

Let us now discuss the specific system studied in this paper, 
i.e., a pair of coupled one-dimensional condensates of inter- 
acting atoms with Hamiltonian H - H\ + H2 + V, where 



Hi 



dx 



Si n 



(1) 



For / = 1,2, n,(.'i:) is the conjugate momentum to bosonic 
field <l>,(x), and the coupling term has a sine-Gordon form 
V = -2^ Jdxcos [(Di(x) - <^2{x)] H. Physically, <l),(jc) and 
n,(x) respectively represent the phase and the density fluc- 
tuations (with respect to a constant background density) of 
condensate ; at position x, and, v, and gi are respectively the 
sound velocity and the Luttinger parameter As seen in Fig.lT] 
- is an effective tunneling amplitude per length (a is a mi- 
croscopic length scale), which can be tuned by changing the 
height of the optical potential barrier. 

A comment on the dimensions of the quantities above is 
in order. We have set fi to unity, and identified the units 
of time and inverse energy. Representing length and energy 
by { and e respectively, the field <l>(x) is dimensionless, its 
conjugate momentum n(x) has dimension and, v and A 
respectively have dimension is and s. Let us now use the 
harmonic approximation (i.e., expand the cosine term around 
^\ix) - ^2ix) - and keep the leading quadratic term). This 
approximation is justified (at least for the initial equilibrium 
state of two coupled condensates) in the limit of large Lut- 
tinger parameters where the cosine term is relevant. As we 
will check a posteriori, although the differences <l>i (x) - <l>2(x) 
typically increase by an optimal evolution designed to cool 
one of the condensates, for some range of parameters, one can 
keep them reasonably small during the evolution so that the 
harmonic approximation remains valid. We can then write the 
Hamiltonian in momentum space as a collection of harmonic 
oscillators: 

/ cj>0 L o' J 

+ ^2A(cDf'-cDf2-)2^s^^3^ 

q>0 

where the superscript % (3) indicates the real (imaginary) 
part. Note that and 11;, respectively have dimension £^^~ 
and '^2. We have not included in Hamiltonian (|2]i the q = Q 

term //„ = ^ 2,- ^ (a^' - K)' + "^if- which is re- 

sponsible for changing the particle number A^' of condensate 

; = 1,2 ((tj) is conjugate to A^' and is the background den- 
sity with L representing the system size) lfl6ll . Our goal here 
is to reduce the energy without reducing the particle number. 



and neglecting Hq is a reasonable approximation if A is not 
too large. Note that to prepare number-squeezed states with 
optimal control, we need to work only with a single-mode 
Hamiltonian Hq (l^ . 

For a given protocol A(f), each mode q in the Hamilto- 
nian (|2]i evolves independently. Although the modes do not 
interact in Eq. they all evolve with the same protocol A(f), 
which induces correlations between them. Therefore, we have 
a fundamentally many-mode problem even without (sublead- 
ing) mode-coupling terms, which we have neglected. The first 
step, however, is to analyze the dynamics of a single mode q 
consisting of just two coupled harmonic oscillators, namely, 
H - Hi + H2 + H\2, where = j(pj/mi + kixf) ^ w/a] a,- and 
H12 - 4(xi - X2)^ with 



— , ki = -Vigiq , 

nVi 7T 



^ _ 4A 
a 



(3) 



Let us assume the initial thermal state is prepared at A - Aq. 
We then evolve the system with a time-dependent protocol 
A(t) (with the constraint < A{t) < /l^ax - 4A„,ax/fl)- For any 
A, we can write the single-mode Hamiltonian as // = jP^P + 



^-X^K{A)X, where = (4i=, ^), 



2 

and 



( yjmixi, ^^m2X2) 



KiA)^ 



(ki+A)/mi —A/ ^jm\m2 
-A/ yOTim2 (k2 + A)/m2 



We can then diagonaUze the above symmetric matrix as 
K{A) - Q{A)Q.{A)Q'''{A), where Q{A) is an orthonormal matrix 
of eigenvectors and Q - diag(dZ<j , ajf), with d), a normal-mode 
frequency. In terms of wi,2(/Io)j the initial density matrix is 
given by 



Po 



— J_g-/Soii'i(-io)oI('lo)fli(-lo)g-;3oM2(/lo)fl2('lo)fl2Uo) 

-2^ 



where aj(Ao) is the annihilation operator for normal-mode 
j = 1,2, which can be written in terms of the annihilation 
operators a, of oscillators ; = 1 , 2 as 



aj = ^J]Qkj{Tjkak+gjkal), 




(4) 
(5) 



For a system evolving with A(t), we can write the Heisen- 
berg annihilation operator of oscillator 1 (or 2) in terms of the 
initial normal-mode operators as 

ai(0 = ^ [ui(t)ai(Ao) + v,(f)a; (io)] , 



where and v, are some complex coefficients satisfying the 
following constraint: 



|Mll' + l«2l'-|vil'-|V2l' = l. 



(6) 
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The initial condition (at f = ) for the coefficients 
above are obtained by inverting Eq. (|4|i, i.e., aj - 

I Tjk Qjk {Tkj cik - Qkj ai). 

To compute m,(t) and v,(t), it is convenient to con- 
sider a piece-wise constant protocol determined by a 
sequence (/I,, f,), for ; - \...n, so that fli(T) - 



mutation relation 



g-iH{Aoti^ Using the com- 



[ai,H{A)] = —(a; - a\) + - } J^ iaj + a)). 



with KiA) = Q'^{Ao)K{A)Q(Ao), we then find that these co- 
efficients at r = r are obtained by integrating the following 
equations of motion from t = to t - t: 



Uj - — 



Kjk{T-t) 



k S^j^k 
> {Uk - Vk) — 1=^ 



+ {uj + Vj)a)j 
- {iij + Vj)a)j 



(7) 
(8) 



where all normal-mode frequencies ca are calculated ?X A - 
Aq. Notice that the equations above depend on the final time 
T, and should not be used for computing m, and v, at f t. 
When working at fixed t, it is helpful to define A{t) = A{t - 1), 
which makes the equations local in time. Finding the optimal 
A immediately yields the optimal A. 

The equations of motion (|7]l and dH), together with their 
initial conditions, uniquely determine a\{T) as a functional 
of A{t). (Notice that the same equations with different ini- 
tial conditions can be used to find 02(1") as well.) Our goal 
is to minimize an appropriate cost function, such as the ex- 
cess energy of oscillator 1, over all permissible controls A{t). 
For a single oscillator, the excess energy is proportional to 
the average number of excitations, which can be written as 
(«i(/)) = tr |^fl[(f)ai(0 po]- In terms of the dynamical variables 
M,(f) and y,(f), the above trace simplifies to 



where n,(0) s tr[fl,(^o)^«;(/lo)pol = - \\\ 

Note that exactly the same formulation describes the many- 
mode problem [Eq. (|2|i]. In this case, we have to multiply 
the number of dynamical variables by the number of modes. 
The equations of motion [Eqs. (|7]i and dHJ] still hold for each 
mode, with parameters depending on q as in Eq. (|3]l. Appro- 
priate many-mode cost functions can be constructed from cost 
functions for individual modes. For example, we can simply 
add (n'j (f)) to obtain the total number of excitations {N\ (t)} in 
condensate 1, or weight them by the mode frequency to find 
the total excess energy (£i(f)) in the condensate: 

<M(0) = 2 2 {n\(t)}, (£i(f)) = 2vi ^ q{nlm (10) 

0<q<A 0<q<A 



where the factor of two accounts for real and imaginary com- 
ponents of Hamiltonian (|2|i and A is a momentum cutoff. Ad- 
ditionally, we may also consider (Ci(f)) - l^o<q<h{^\it)) I q, 
which is relevant for enhancing the fringe contrast of matter- 
wave interferometry experiments UtIi . 

The problem we have formulated thus far is a typical prob- 
lem in optimal control theory applied to quantum dynam- 
ics 1 1 8I - I22II : we have a set of dynamical variables with 
given initial conditions (m, and v, in our case), which evolve 
with given equations of motion [Eqs. (|7|i and ([8])] that depend 
on some admissible control parameter(s) (0 < A{t) < /Imax)- 
The challenge is to find an admissible optimal control such 
that a given cost function of the dynamical variables [Eq. ( fTOl ) 
in our case] is minimized at a given time t. 

Let us now turn to the main questions of this work: What 
do the optimal A{t) protocols look like? How can we find 
them? How much can they cool a system? Using Pontryagin's 
maximum principle, we argue that optimal protocols are bang- 
bang, i.e., A{t) is either zero or equal to /Imax at any given time. 
As mentioned earlier, we demonstrate that a direct simulated- 
annealing calculation can yield these optimal protocols. We 
also find that, depending on the parameters of the problem, it 
is possible to significantly cool down one of the condensates. 

Let us now briefly review Pontryagin's maximum princi- 
ple. Consider a set of dynamical variables {x(f)) that satisfy 
the equations of motion Xj - fji{x, a}), with Xj{Q) - x^j, for 
a set of admissible controls {a{t)}. The goal is to maximize a 
payoff function g{{x{T)}) over all such {a(t)}. The key to Pon- 
tryagin's maximum principle is the following optimal-control 
Hamiltonian: 

mix, p, a}) = Yj Pj(f^ //(^' '^D' (11) 

where pj{t) is a "momentum" conjugate to Xj{t). The Pon- 
tryagin's theorem states that for the optimal control {a*{t)}, 
and the corresponding {x*, p*}, we have 

Jif* = J^{{x\p\a}) = maxJ^({x*,p\a}), (12) 



with bound- 

dx. 



where x and p satisfy i* - and p* - 

ary conditions x*(0) - x^j and /'*(t) = j^.g({x*(T)]). It is now 
easy to observe that since, for all modes q, Eqs. (|7]) and ([8]) 
are linear in A{t), the optimal-control Hamiltonian [Eq. ( fTTT il 
is also a linear function of A{t) in our case. We then immedi- 
ately deduce from Eq. ( fT2] i that, unless Jif* identically van- 
ishes over a finite time interval, the control A{t) can only take 
two values, namely, zero and /Imax- 

The Pontryagin equations are not easy to solve numerically 
for many modes. We thus use our direct MC method, with- 
out utilizing any assumptions regarding the bang-bang nature 
of the protocol. The kinetics of the simulated annealing con- 
sist of varying a randomly chosen /I, (of the discretized pro- 
tocol) by a small random amount. As found in Ref. 11411 . 
such simulations converge very well in the number of dis- 
cretization points. To compute the cost function (in each 
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MC step), it is convenient to define new dynamical variables 
= ^JLbj{AQ) (uj + Vj^ and Oj = (uj - vy) / ^I^JJIq), which 
satisfy \^) - -iK{A)\6) and \0) - -i\(p) in matrix notation (this 
change of variables allows us to diagonalize 2x2 matrices 
instead of 4 x 4 ones). By solving the above equations in 
terms of the eigenvalues and eigenvectors of the 2x2 matrix 
K(Am) - K{A„-m+i), we can then write simple recursion rela- 
tion for 10) and \9), which yield their values at time r after n 
iterations. 

One comment is in order before proceeding. In addition 
to optimizing over A(f), we have the freedom to choose the 
initial A^ anywhere between zero and /l,nax- There is a rigor- 
ous lower bound on {iii) [Eq. (|9]l], which follows from con- 
straint Q: <«i(f)> > Hmin = min««i(0)), <«2(0))). We can 
show that «niin is a decreasing function of Aq. This suggests 
that it may be advantageous to set Aq - /Imax- Although the 
actual cost functions we are able to reach by our minimization 
procedure are typically much larger than this lower bound, by 
trying several values of Aq in our numerics, we have found that 
the best cooling is in fact achieved for Aq - /l,nax- In Fig.|2^, 
we show a typical protocol obtained by MC simulations. We 
converge to a bang-bang protocol by an unbiased simulation, 
which samples all the intermediate values of A, and, a pri- 
ori, does not assume anything about the shape of the protocol. 
Surprisingly, minimizing N\, &i, or Ci leads to very similar, 
albeit nonidentical, protocols (here we show the protocol ob- 
tained by minimizing Ni). To further check the consistency 
with the Pontryagin's theorem, we also computed the deriva- 
tive (9^ J^, the sign of which determines A{t) through Eq. ( fT2] i. 
for cost function {Ni{t)). 

To construct J^, we need to treat the real and imaginary 
parts of <pj and 6j as separate dynamical variables with their 
own conjugate momenta. We can then construct a complex 
variable tt^, whose real (imaginary) part is the conjugate mo- 
mentum to the real (imaginary) part of (pj, and similarly for 9j. 
For each q, we then have I.tt'^) = -i\n'') and I.tt") = -;A'(3)|7r'*>. 
The boundary conditions at f = t depend on the cost func- 
tion [see the boundary conditions below Eq. (fT2] i1 and, for 
<M(t)), can be written as 7rf(T) = 0,(t) - (2<«,(0)) + D sjo) 
and 7rf(T) = 0,(t) - (2<n,(0)> + l)w,(0) 0i{T). Given a pro- 
tocol A{t), we can solve for (p and 6 forward in time, con- 
struct 7r'*(T) and n^(T) from the boundary conditions above, 
solve for tt*^ and tt" backward in time and finally construct 
= Yjq{^'q\^AKq{A)\Oq), which immediately yields d,\-^- 
The results are shown in Fig. |2^, and show excellent agree- 
ment with the simulations. 

In Fig. lib, we show how the cost function {N\{t)) changes 
when evolving with the optimal protocol. An interesting fea- 
ture of the evolution is that '^^^'^^'^^ just before quench- 
ing to A{t) - 0. Keeping the subsystems coupled would do 
a better job in reducing {N\ (t)) locally (in time) but would 
not lead to global optimization in total time t. We can also 
check the harmonic approximation a posteriori by computing 
J J dx{[<S>i{x) - (S>2(x)]~). We find that as long as the approx- 
imation is valid initially, and /t^ax is large enough, this quan- 



L/a = 32,jSoT = 1.5, Amax^ = 5 




FIG. 2: a) Typical protocols obtained with unbiased simulated an- 
nealing, and the derivative of the optimal control Hamiltonian with 
respect to control A, whose sign determines the protocol. The simu- 
lations converge to bang-bang protocols in excellent agreement with 
sgn b) Reduction of <M(f)> due to evolution with the two 

optimal protocols above. Cooling condensate 1, may also cool down 
the other condensate for free. For large Luttinger parameters, the ar- 
gument of the cosine term remains much smaller than one during the 
evolution, i.e., the harmonic approximation remains valid. 



tity remains smaller than one and the harmonic approximation 
holds throughout the evolution (if the system is not too long, 
the spatial fluctuations of [Oi(x) - <S>2{x)]^ are small). Inter- 
estingly, the optimal protocol designed for reducing the en- 
ergy of condensates 1 turns out to also cool down condensate 
2. This is not a violation of the second law of thermodynam- 
ics as our process is not cyclic: we start from two coupled 
condensates with H - H\ + H2 + V , and end up with two de- 
coupled ones with H - Hi + H2- The process only reduces 
the expectation value {Hi + H2), while {Hi + H2 + V), which 
corresponds to the initial Hamiltonian, actually increases. 

The eff'ective cooling described here is an out-of- 
equilibrium reduction of the excess energy, and does not im- 
ply thermal equilibrium. If the low-energy system equilibrates 
afterwards, however, it will have a lower temperature. To di- 
rectly bring the system close to thermal equilibrium, one can 
instead minimize the trace distance between the density ma- 
trix and thermal density matrices at varying target tempera- 
tures, and find a balance between a small trace distance and 
a low target temperature. We do not pursue this approach 
here because, under realistic circumstances, each decoupled 
condensate is expected to eventually decohere, and reach an 
effective temperature determined by its excess energy 1231]. 

Finally, we discuss the performance of our optimal proto- 
cols, which depends on several (diemnsionless) parameters, 
including the two Luttinger parameters, /Jqt, v/t/a, and Lja. 
However, Amax"!" and the ratio V2/V1 seem to have the most 
pronounced effect on the performance (measured by the ratio 
of the achieved energy to equilibrium energy at Pq). Typi- 
cally, the energy can be reduced by a factor of 3 to 5 when 
V2/V1 is of order unity. For a highly asymmetric system with 
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V2/vi = 100, we achieved an energy reduction by a factor of 
40 with a system size of L/a = 64 and other dimensionless 
parameters of order unity. 

Let us now comment on a possible extension of our scheme 
to arbitrary systems. To perform our MC simulations, we need 
to be able to efficiently compute desired cost functions for any 
allowed protocol (which is the case for our system in the har- 
monic approximation). Cooling down more complicated sys- 
tems such as the two-dimensional fermionic Hubbard model 
(or even our system in regimes where the full sine-Gordon 
term is needed) is of considerable interest for quantum sim- 
ulations. The generality of our MC method, however, raises 
an intriguing possibility for a universal approach: if one can 
automate the processes of system initialization (initial cool- 
ing), unitary evolution (with a tailored protocol), and mea- 
surement of the figure of merit (e.g., energy), then the system 
itself can be used to perform such MC simulations. The cost 
function can be measured (instead of computed), and then fed 
into the MC algorithm. This would provide a powerful means 
of preparing desired states in arbitrary systems, and may open 
the door to the quantum simulation of unsolved condensed- 
matter models. Such integration of experiment and simula- 
tion has in fact been applied to the control of chemical reac- 
tions II24II . and more recently to some aspects of cold atom 
experiments 1125 1 . 

In summary, we demonstrated that nonadiabatic optimal 
control of quantum evolution can be used to push the bound- 
aries of atomic cooling. Contrary to the conventional asso- 
ciation of nonadiabatic effects with heating, we showed that 
breaking away from the adiabatic limit, in a controlled way, 
can in fact help cool down quantum systems. We applied 
this idea to a system of two coupled elongated condensates. 
Through simple and direct MC simulations, we found optimal 
protocols which agree with theoretical predictions based on 
Pontryagin's maximum principle, and are effective in reduc- 
ing the excess energy. Such MC simulations can be potentially 
performed by the system itself giving access to a universal 
cooling scheme. 
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